subroutine ellipse(MatrixS,MatrixA,OldOval,NewOval,T)
implicit none

type :: Oval
  real :: l(2,1),r(2,1)
  real :: a2,b2
end type Oval

real :: MatrixS(2,2), matrixA(2,2)
type(Oval) :: OldOval,NewOval,Coord1,Coord2
real :: C(2,2),L(2,2),M(2,2),D(2,2),Center(2,2),Tmp1(2,2),&
& Tmp2(2,2),Tmp3(2,2)
real :: vec_i(2,1),vec_j(2,1)
real :: T
real :: omega


vec_i(1,1)=1
vec_i(2,1)=0
vec_j(1,1)=0
vec_j(2,1)=1
!--------------------------------------------
Center(1,1)=OldOval%a2
Center(1,2)=0
Center(2,1)=0
Center(2,2)=OldOval%b2
!----------------------------------------------
call eig(MatrixS,Coord1%a2,Coord1%b2,Coord1%l,Coord1%r)
L(1,1)=1+Coord1%a2*T
L(1,2)=0
L(2,1)=0
L(2,2)=1+Coord1%b2*T
call turn(OldOval%l,Coord1%l,C)

!-----------------------
omega=2*MatrixA(1,2)
M(1,1)=Cos(-omega*T)
M(1,2)=Sin(-omega*T)
M(2,1)=-Sin(-omega*T)
M(2,2)=Cos(-omega*T)
!-------------------
Tmp1=matmul(L,C)
Tmp1=matmul(M,Tmp1)
call inv(Tmp1,Tmp2)
Tmp3=matmul(Center,Tmp2)
Tmp3=matmul(transpose(Tmp2),Tmp3)
!---------
call eig(Tmp3,Coord2%a2,Coord2%b2,Coord2%l,Coord2%r)

!-------------- 
call turn(Coord1%l,vec_i,D)
NewOval%l=matmul(D,Coord2%l)
NewOval%r=matmul(D,Coord2%r)
NewOval%a2=Coord2%a2
NewOval%b2=Coord2%b2
end subroutine ellipse
!==================================================================================================================================
subroutine inv(matrixM,result_m)
 real :: matrixM(2,2),adjM(2,2)
 real :: result_m(2,2)
 real det
 !call printM(matrixM)
 det=matrixM(1,1)*matrixM(2,2)-matrixM(1,2)*matrixM(2,1)
 !write(*,"('DET in Function=',E20.8)") det
 !!if det .EQ. 0
 !!!
 adjM(1,1)=matrixM(2,2);adjM(1,2)=-matrixM(1,2);adjM(2,1)=-matrixM(2,1);adjM(2,2)=matrixM(1,1)
 !call printM(adjM)
 result_m(1,1)=adjM(1,1)/det;result_m(1,2)=adjM(1,2)/det;result_m(2,1)=adjM(2,1)/det;result_m(2,2)=adjM(2,2)/det
return 
end
!==================================================================================================================================
function det(matrixM)
real :: matrixM(2,2)
real det
det=matrixM(1,1)*matrixM(2,2)-matrixM(1,2)*matrixM(2,1)
end
!==================================================================================================================================
subroutine printM(matrixM)
real :: matrixM(2,2)
write(unit=*,FMT="('|',E20.8,E20.8,'|'/'|'E20.8,E20.8,'|')") matrixM(1,1),matrixM(1,2),matrixM(2,1),matrixM(2,2)
end
!==================================================================================================================================
subroutine normalize(vec,nvec)
real :: vec(2,1),nvec(2,1)
real norm
norm=Sqrt(vec(1,1)**2+vec(2,1)**2)
nvec(1,1)=vec(1,1)/norm
nvec(2,1)=vec(2,1)/norm
end
!==================================================================================================================================
subroutine eig(M,lamda1,lamda2,vec1,vec2)
real :: M(2,2)
real :: vec1(2,1),vec2(2,1),v1(2,1),v2(2,1)
real lamda1,lamda2,Delta
real det

!if det.LE. else
!
!

Delta=(M(1,1)+M(2,2))**2-4*det(M)
lamda1=((M(1,1)+M(2,2))+Sqrt(Delta))/2
lamda2=((M(1,1)+M(2,2))-Sqrt(Delta))/2
If (Delta .GE. 1E-7) then
call value2vec(lamda1,M,v1)
call value2vec(lamda2,M,v2)

call normalize(v1,vec1)
call normalize(v2,vec2)
Else
vec1(1,1)=1
vec1(2,1)=0
vec2(1,1)=0
vec2(2,1)=1
Endif

Contains
	subroutine value2vec(val,MatrixTMP,vec)
	real val
	real :: MatrixTMP(2,2)
	real :: vec(2,1)
	real val1,val2
	val1=Abs(MatrixTMP(1,1)-val)+Abs(MatrixTMP(1,2))
	val2=Abs(MatrixTMP(2,1))+Abs(MatrixTMP(2,2)-val)
	if (val1.LT.val2)then
		vec(1,1)=val-MatrixTMP(2,2);vec(2,1)=MatrixTMP(2,1)
	else!val1>=val2
		vec(1,1)=MatrixTMP(1,2);vec(2,1)=val-MatrixTMP(1,1)
	end if
	end subroutine

end
!==================================================================================================================================
subroutine cross_vec(a,b,outarg)
real :: a(2,1),b(2,1),m(2,2)
real det
m(1,1)=a(1,1);m(1,2)=a(2,1);m(2,1)=b(1,1);m(2,2)=b(2,1)
outarg=det(m)
end
!==================================================================================================================================
subroutine turn(vec_old,vec_new,M)
real :: vec_old(2,1),vec_new(2,1),M(2,2)
real c,s
!c=dot_product(vec_old,vec_new)
!s=cross_product(vec_old,vec_new)
call dot_vec(vec_old,vec_new,c)
call cross_vec(vec_old,vec_new,s)
M(1,1)=c;M(1,2)=s;M(2,1)=-s;M(2,2)=c
end
!==================================================================================================================================
subroutine dot_vec(a,b,outarg)
real :: a(2,1),b(2,1)
real outarg
outarg=a(1,1)*b(1,1)+a(2,1)*b(2,1)
end

